A simple model of COVID-19 explains disease severity and the effect of treatments

Considerable effort has been made to better understand why some people suffer from severe COVID-19 while others remain asymptomatic. This has led to important clinical findings; people with severe COVID-19 generally experience persistently high levels of inflammation, slower viral load decay, display a dysregulated type-I interferon response, have less active natural killer cells and increased levels of neutrophil extracellular traps. How these findings are connected to the pathogenesis of COVID-19 remains unclear. We propose a mathematical model that sheds light on this issue by focusing on cells that trigger inflammation through molecular patterns: infected cells carrying pathogen-associated molecular patterns (PAMPs) and damaged cells producing damage-associated molecular patterns (DAMPs). The former signals the presence of pathogens while the latter signals danger such as hypoxia or lack of nutrients. Analyses show that SARS-CoV-2 infections can lead to a self-perpetuating feedback loop between DAMP expressing cells and inflammation, identifying the inability to quickly clear PAMPs and DAMPs as the main contributor to hyperinflammation. The model explains clinical findings and reveal conditions that can increase the likelihood of desired clinical outcome from treatment administration. In particular, the analysis suggest that antivirals need to be administered early during infection to have an impact on disease severity. The simplicity of the model and its high level of consistency with clinical findings motivate its use for the formulation of new treatment strategies.

www.nature.com/scientificreports/ simple and yet adequately captures the main clinical findings. Many within-host models of SARS-CoV-2 infection have been published [11][12][13][14][15][16][17] with a wide range of model complexities (from less than 10 to more than 80 parameters). These models mostly differ in terms of the complex interactions between the viral infection and the immune response that they include in their formulation (see Discussion for details). The large difference between model formulations suggests that identifying the key elements having an impact on clinical outcomes is a difficult task. The model we formulated focuses on cells that trigger inflammation through molecular patterns: infected cells carrying pathogen-associated molecular patterns (PAMPs) and damaged cells producing damage-associated molecular patterns (DAMPs). We show that the clearance rate of infected and damaged cells by the innate immune response is of the utmost importance to reach a state of resolved inflammation. Our model can explain the following findings: (i) severe COVID-19 tends to be accompanied by hyperinflammation, (ii) those with severe COVID-19 generally experience a similar viral trajectory as mild cases, albeit with a slower viral load decay after the peak, (iii) the complex and conditional effect of antivirals and corticosteroids on disease severity, (iv) an inefficient type-I IFN response is associated with severe COVID-19, and (v) generation of bystander cell damage and infective removal of these cells are a critical component of severity. Note that this last point is reminiscent of clinical observations that, for example, less cytotoxic NK cells and higher levels of neutrophil extracellular traps (NETs) are associated with severe COVID-19. Overall, the simplicity of the model we propose along with its high level of consistency with clinical observations suggest it is an adequate framework for the study of COVID-19 pathogenesis and the effect of therapy.

Methods
The model. Our goal was to formulate a model that is simple enough to guide intuition, yet complex enough to allow relating with clinical outcomes. A schematic representation is provided in Fig. 1. The model is described below.
(1) dT dt = −βTV − κ 0 TI + δ R R at rate κ 0 I (we assume target cells are exposed to a concentration of type-I IFN that is proportional to I and that puts the cells into an antiviral state). Resting innate immune cells ( D 0 ) become activated ( D 1 ) at rate σ (I + J), where I and J are the number of infected and damaged bystander cells, respectively (we assume the extent of PAMP and DAMP signaling is proportional to I+J) . Also, activated immune cells ( D 1 ) die at per capita rate δ D . Damaged bystander cells ( J ) are generated from an extensive proinflammatory response at a rate that is a Hill function of the number of activated immune cells D 1 . Infected cells die due to viral cytopathic effects at rate δ I and damaged cells die from their injury at rate δ J . The clearance of these cells also occurs by the action of activated innate immune cells at rate κ 1 D 1 . The effect of the adaptive immune response is modeled by adding a constant term κ 2  . We assumed this process follows mass action with rate βV . Infected cells die at constant rate δ I due to viral cytopathic effects. Infected cells produce virions at rate p , which are cleared at per capita rate c.
Modeling inflammation and the innate immune response. The model focuses on molecular patterns that initiate an inflammatory response. Pathogen associated molecular patterns (PAMPs) trigger an inflammatory response by signaling the presence of pathogens 18 . An inflammatory response can also be initiated in the absence of pathogenic infection through recognition of damage associated molecular patterns (DAMPs) 19 . Various conditions promote DAMP expression including hypoxia, low levels of glucose and amino acids, exposure to heat, physical stress or exposure to toxic molecular products [19][20][21][22] . Many of the same conditions that promote DAMP expression can be induced by inflammation itself, for example, an excessive presence of neutrophil extracellular traps (NETs) released during inflammation is linked to immunothrombosis 23 . This can contribute to hypoxia and nutrient deprivation that can lead to DAMP expression 24 . Overall, a positive feedback loop between inflammation and DAMP expression may be triggered 22,24 . Interestingly, patients with severe COVID-19 have increased levels of NETs 25,26 . NETs and platelet dysregulation may be specific to COVID-19 and both are associated with lung microthrombi, conditions that can ultimately favor an inflammatory response in the lungs 27 . The inflammatory response signals a need to eliminate infected or damaged cells and protect the rest of the organism from the perceived danger. Infected cells carrying PAMPs and uninfected cells producing DAMPs are represented by I and J in the model, respectively. In the organism, these molecular patterns promote the downstream recruitment and activation of numerous innate immune cells such as neutrophils, monocytes, macrophages and natural killer (NK) cells 28 . For simplicity, we lumped all cells that take part directly or indirectly (through cytokine signaling) in cytolytic or phagocytic activities into two model compartments, D 0 and D 1 representing resting and activated phenotypes, respectively. In our model, the rate at which D 0 cells become D 1 is proportional to the amount of PAMPs and DAMPs, themselves assumed proportional to the number of cells carrying these molecules, i.e. σ (I + J). We further assumed resting cells D 0 have a homeostatic level D 00 , that they maintain through the recruitment of resting immune cells over the course of infection at rate . We assumed D 0 cells are long-lived compared to the duration of acute infection, and that activated cells decay at rate δ D . The role of D 1 in our model is to promote cytolytic and phagocytic activities to rid the system of PAMPs and DAMPs. Due to the innate nature of the modeled response, we assumed the effect to be similar on the decay of both I and J at per capita rate κ 1 D 1 . This was partly motivated by the behavior of NK cells, which target injured epithelial cells through stimulation of receptor NKG2D, as well as infected cells having down-regulated major histocompatibility complex class 1 (MHC-I) molecules or upregulated MHC class I polypeptide-related sequence A (MICA) or sequence B (MICB) 29,30 .
We assumed that high levels of inflammation can induce significant cellular stress and promote the expression of DAMPs in bystander cells. This is represented by the Hill term in Eq. (5) for the rate of generation of uninfected cells expressing DAMPs with ν being the maximum damage rate, K being the D 1 concentration leading to half of the maximum damage rate, and m is the Hill coefficient. We assumed damaged cells decay at rate δ J < δ I in the absence of an innate immune response.
Finally, we assumed an additional innate immune response, e.g., a type-I interferon response. Type-I interferons signal danger to neighboring cells, making the latter refractory to infection ( R) 31 . For simplicity, we approximated the amount of interferon in the microenvironment to be proportional to the number of infected cells 32 . Accordingly, we modeled the rate at which susceptible cells become refractory (R) by κ 0 I. Refractory cells expected to naturally revert to a susceptible state after 1/δ R .
Modeling the adaptive immune response. The model mainly focuses on the innate response. To avoid unnecessary complexity, the adaptive immune response is represented by a single term that includes two parameters: κ 2 , www.nature.com/scientificreports/ which is the maximum decay rate of infected cells due to the adaptive immune response, and τ, which represents the time post-infection that the adaptive response takes effect. This approach is similar to that in Pawelek et al. 33 . This response is thus modeled using the indicator function I t>τ (t).
Structural analysis of the model. A substantial minority of individuals with COVID-19 experience a state of sustained high inflammation or hyperinflammation 9 . We evaluated whether the model allowed the existence of such a hyperinflammatory state. We used D 1 as a marker of inflammation. Through a bifurcation analysis, we searched for the existence of a stable steady state of resolved inflammation ( D 1 = 0) and the existence of a second stable steady state of hyperinflammation where D 1 is maintained well above 0. This type of analysis also allows identifying basins of attraction, i.e. regions in the space of variables that lead to specific inflammation trajectories over time. We used analytical techniques combined with the numerical bifurcation software Matcont 34 , a Matlab software package designed for the analysis of equations such as Eqs.
Set of virtual markers. To analyze model behavior, a set of virtual markers were computed. These markers were used to compare model predictions to clinical observations and to investigate determinants of severe COVID-19. These markers are: (i) the peak viral load (peak VL, maximum value of V); (ii) the time of peak viral load; (iii) the difference between the viral load at its peak and the viral load 5 days after the peak; (iv) peak D 1 (used as a proxy measure for peak inflammation); (v) the time of peak D 1 ; (vi) D 1 at 60 days post-infection. We further computed the hyperinflammation index, which takes a value of 1 if the value of log 10 (D 1 ) at the end of the simulated period was more than 99% of its peak value, and 0 otherwise. This resulted in the categorization of each simulation into inflammation trajectory groups: (i) resolved inflammation (R) or (ii) hyperinflammation (H). Finally, we computed the Disease Score, defined as the total number of cells (I + J) that died over the 60-day period post-infection. This score is meant to describe disease severity, with higher scores representing more severe COVID-19.
Identifying a space of realistic parameter values. The next step consisted in identifying a space of parameter values for which model predictions are consistent with a minimal set of clinical observations to allow in silico investigation of the model. In particular, acceptable parameter values should result in: (i) peak viral loads (VL) values between 4 and 10 log 10 , (ii) peak viral loads achieved 2 to 14 days following infection, and (iii) an innate immune response during the course of infection (the activation of at least 1% of all resting innate immune cells). These are referred to as conditions (i)-(iii). Condition (i) was chosen to represent peak viral loads observed using nasopharyngeal swabs in a population of individuals infected by the virus 7,37 . Condition (ii) was chosen based on reports of peak viral loads occurring around the time of symptom onset and symptom onset primarily occurring within 14 days of infection (median 4-5 days) 14,16,32,38 . Finally, we chose condition (iii) to ensure D 1 reaches high enough values for an observable effect of D 1 on either I(t) or J(t).
The determination of the space of parameter values that satisfy the above conditions was done iteratively. First, we established bounds of values for each parameter using literature estimates when available, results from the bifurcation analysis to ensure that hyperinflammation was achievable, and preliminary simulations when no estimate was available (see Table 1). We then selected n = 100,000 vectors of parameters from the parameter space using a Latin hypercube approach to minimize the chances of unexplored multidimensional subspaces 39 . Simulations were performed for each selected vector of parameters using initial conditions listed in Table 2 to predict infection and inflammation trajectories from day 0 (day of infection) to day 60. We identified regions in the parameter subspace leading to unacceptable results based on the conditions (i)-(iii) listed above. We subsequently refined the parameter space and repeated the procedure until we were satisfied that a randomly selected vector of parameters would lead to a high likelihood of satisfying acceptance criteria (i) to (iii) (> 70% acceptance probability). Table 1 describes all parameters, the resulting space of parameter values along with references from the literature when applicable.
Simulating clinical observations. We sampled a large number of sets of parameter values (n = 1,000,000) from the space of parameter values described in Table 1 using a Latin hypercube approach 39 . Accordingly, the marginal distribution of values being sampled for each parameter was uniform across the considered range (Table 1). For each selected vector, simulations were performed using initial conditions listed in Table 2 to predict infection and inflammation trajectories from day 0 (day of infection) to day 60. Analyses were only performed on the simulations satisfying the acceptance criteria (conditions (i)-(iii)). We first investigated the distribution of parameter values and virtual markers using histograms. We used violin plots to study bivariate associations between parameters and inflammation trajectory groups (resolved inflammation or hyperinflammation). Multivariate analyses were performed using decision trees 45 . We used the hyperinflammation index as the variable being predicted by model parameters. We first obtained a single tree. Cross-validation and deviance plots were used to guide the choice of the optimal tree 45 . To account for uncertainty around the formulation of a single tree, we also performed a Random Forest analysis, generating 100 trees, and reported the mean decrease in GINI index (a measure of decreased node impurity from choosing a parameter for tree splits, i.e. its ability to discriminate inflammation trajectory groups) 46 .
Simulating the effect of treatments on COVID-19 severity. Finally, we evaluated if the model could replicate clinical findings regarding the treatment of COVID-19. For this purpose, we selected 10,000 simulations leading to lower Disease Severity scores and the same number of simulations leading to higher scores from the sample of accepted simulations. To allow comparison with clinical data, we assumed the former group www.nature.com/scientificreports/ represents mild/moderate disease, while the latter represents severe COVID-19. To simulate the use of corticosteroids, we assumed a reduction of D 1 by 50% for a period of 10 days. To simulate the use of potent antivirals, we decreased parameters β or p to 1% of their original value at peak viral load for the remainder of the simulation. For each simulation, we modified parameter values at the time of peak viral load as it is estimated that the peak is reached within the few days following symptom onset 14,32,38 . We also simulated treatment administration a day prior to the expected viral load peak to study the effect of early treatment. The difference in log 10 Disease Score (with vs without treatment) was computed and reported.  40 . A larger range of values was used to account for variability between individuals and uncertainty in parameter values due to a possible correlation with parameter β . Accordingly, it was inefficient to sample β and p independently. It was determined through linear regression that a value of log 10 p determined by 6.8-1.0*log 10 β + U, where U is a uniform random variable with range [ 0.4;0.5], greatly enhanced the likelihood of peak viral loads being between 4 and 10 log 10 and occurring between day 2 and 14 post-infection . A value of the order of 1 × 10 -6 was used in Jenner et al. 11 . The range of values was adjusted to ensure that at least 1% of all resting innate immune cells activated over the course of infection and that high inflammation levels ( D 1 > 1 × 10 6 ) did not occur prior to peak VL . Interestingly, the analysis shows that under certain biologically plausible parameter values, there exists bistability in the system with both a stable hyperinflammatory state and a stable resolved inflammatory state. In more technical terms, there are three equilibria in the model: a stable high-inflammation state corresponding to hyperinflammation, an unstable equilibrium with non-zero, but low, inflammation and another stable equilibrium corresponding to resolved infection and resolved inflammation. The first two steady states appear/disappear following a saddle-node bifurcation as shown in the bifurcation plots (Figs. S1 and S2). The bifurcation analysis identified important parameters that dictate the existence of a hyperinflammatory steady state: (i) κ 1 , which represents the effect of the innate response ( D 1 ) on the clearance of cells carrying PAMPs and DAMPs, and (ii) the parameter ν that dictates the amount of bystander cell damage resulting from inflammation (Fig. S1). Our analysis shows there is a threshold value of κ 1 beyond which the hyperinflammatory state ceases to exist (Fig. S1A). Further, hyperinflammation can only occur if there is enough immune driven inflammation beyond a threshold value of ν (Figs. S1B and S2B).
When bistability exists in the system, inflammation trajectories could either converge to one stable state or the other over time. What determines the long-term inflammation trajectory is the amount of bystander cell damage J caused by inflammation over the course of the infection. Mathematically, this is represented by a saddle-node bifurcation. The unstable lower branch of equilibria acts as a separatrix between the hyperinflammatory and resolved-inflammation states and so defines the basin of attraction of the hyperinflammatory state (see Fig. S1). Thus, hyperinflammation can only occur if the infection induced inflammation is severe enough to force the system across this separatrix and into this basin.
Overall, these results suggest that the clearance of cells carrying PAMPs and DAMPs and the amount of bystander cell damage from inflammation are important determinants of disease outcomes. When the clearance rate of infected and damaged cells is sufficiently high or bystander cell damage is low, infection leads to non-severe outcomes.

Model simulations of a virtual cohort of infected individuals.
We next simulated the model by sampling n = 1,000,000 sets of parameter values across biologically plausible ranges (Methods). We use these simulation results to analyze the different viral load and inflammatory response trajectories in the population, such that key determinants of disease outcomes can be further identified. In our analysis, we only included simulations satisfying the acceptance criteria that are consistent with broad patterns seen in clinical studies (described in the Method section as conditions i-iii) (accepted simulations, n = 739,465 of 1,000,000 or 73.9%). The distribution of parameter values that led to accepted simulations are presented in Supplementary Material (Fig. S3). Most parameters preserved their sampling distributions, i.e., the distribution remained uniform across the considered ranges after discarding simulations that do not satisfy the acceptance criteria. Further, there was little correlation between almost all parameters. However, there was a strong correlation between β and p ( 0.99). Consequently, parameters β and p could not be independently sampled in order to output acceptable simulations (Fig. S4). Table 1 describes the sampling strategy used to account for this dependency.
A strong association between hyperinflammation and the disease score. To represent disease severity, we define the Disease score as the total number of infected and bystander cells I and J that died over a 60-day period. Simulations that give rise to the Disease Score distribution are illustrated in Fig. 2. They are categorized in groups by inflammation trajectory: (i) the inflammation marker D 1 decreased after peak inflammation or (ii) the inflammation marker D 1 kept increasing or was maintained at high level following peak inflammation (hyperinflammation). There was a direct link between the presence of hyperinflammation and the Disease Score (see Fig. 2a). Overall, 13.5% of accepted simulations exhibited hyperinflammation. Hyperinflammation and a high Disease Score were both associated with a higher number of cells that died following injury from inflammation (Fig. 2b). www.nature.com/scientificreports/

Similar viral load but divergent inflammation trajectories. Viral load and inflammation trajectories
over time by inflammation groups are presented in Fig. 3. Simulations resulting in both resolved inflammation or hyperinflammation exhibited similar viral load dynamics. The distribution of peak viral load and the time to reach this peak after infection largely overlapped between groups. However, those with resolved inflammation tended to have a faster VL decay after peak VL (Fig. 3c). The most remarkable difference between groups pertained to the dynamics of the inflammation marker D 1 . Peak levels of inflammation were higher for simulations resulting in hyperinflammation. For those with resolved inflammation, peak inflammation was generally observed around the time of peak viral load while for those with hyperinflammation, peak inflammation was observed much later (Fig. 3h).
Association between hyperinflammation and characteristics of the innate immune response. Next, we compared the distribution of parameter values between inflammation trajectory groups. Figure 4 shows violin plots for the 15 parameters that were allowed to vary between simulations. The most striking differences between groups are observed for parameters κ 1 , , κ 0 , σ , δ D and ν . Lower values of κ 0 (type-I IFN secretion and/or response) and κ 1 (cytolytic and phagocytic activities of innate immune cells) were associated with a greater risk of hyperinflammation and a high Disease Score. Similarly, higher values of (innate immune    www.nature.com/scientificreports/ Prediction of the risk of hyperinflammation from characteristics of the innate immune response. We used regression tree analysis to reveal the discriminatory importance of parameters from a multivariate perspective. The regression trees attempted to discriminate simulations leading to hyperinflammation from those leading to a resolved inflammatory state. Results are shown in Fig. 5. The tree has a root (left side), branches, nodes (where branches separate) and terminal leaves (right side). Simulations enter at the root and separate at nodes. At each node, one parameter and one threshold value were chosen by the regression algorithm based on their ability to separate simulations into more homogeneous groups in terms of inflammation trajectory. Simulations are directed toward the lower branch if the value for the chosen parameter for the simulation is smaller than the threshold, or toward the upper branch otherwise. This resulted in 5 terminal leaves that partitioned all simulations into 5 groups. Longer branches following a node signify that the node allowed a better separation of the two inflammation trajectory groups. Figure 5 illustrates the multivariate conditions favoring hyperinflammation. The most important parameters pertain to type-I IFN response, the ability of the system to clear cells carrying PAMPs and DAMPs ( κ 0 and κ 1 ) and the immune cell activation and recruitment rate ( σ and ). In particular, combined low values of κ 0 and κ 1 and high values of σ and led to a dramatically increased chance of hyperinflammation (73.2% risk of hyperinflammation).

Differences in the effect of treatments in terms of treatment type, administration time and predicted disease severity.
To validate and demonstrate the utility of our model, we simulated the use of corticosteroids and antivirals in infected individuals and compared the model results with clinical findings. First, we modeled corticosteroid treatment by assuming the treatment leads to a reduction of D 1 by 50% for a period of 10 days following peak infection. The in silico administration of corticosteroids had a remarkably different effect on the Disease Score depending on the inflammation trajectory groups (Fig. 6). Among those simulations where resolved inflammation is predicted in the absence of treatment (orange plot), corticosteroids were often detrimental (32% chance of an increase in Disease Score). Such a detrimental effect was not observed among those for which hyperinflammation was predicted in absence of treatment (pink plot). Greater improvements were observed in the latter group: 23% had a greater than 0.5 log 10 decrease in Disease Score, compared to only 0.7% in the former group. There were only small differences in the effect of corticosteroids across the investigated timing of drug administration (at peak viral load in Fig. 6a and a day prior in Fig. 6b). Corticosteroid treatment also had an effect on the slope of viral load decay, which was substantial for those simulations of mild/moderate disease and lesser among simulations of severe COVID-19 (Fig. S8).
Simulations of the effect of antivirals were performed by decreasing parameters p or β to 1% of their original value at peak viral load. Modifying p or β had a very similar effect on Disease Scores . None of the simulations revealed a substantial increase in Disease Score. However, the impact of antivirals was very different between inflammation trajectory groups and between times of drug administration. Earlier drug administration led to large improvements in Disease Score for both groups (Fig. 6b). However, administration at peak viral load www.nature.com/scientificreports/ led to much reduced improvements in Disease Score (Fig. 6a), in particular for those simulations leading to hyperinflammation in absence of treatment. In this case, infection had already driven inflammation to high levels by the time peak viral load was reached (although peak inflammation was not reached until much later), which translated into a number of damaged cells that corresponded to being within the basin of attraction of hyperinflammation (Fig. S1).

Discussion
In this work, a mathematical model was formulated to represent the within-host dynamics of COVID-19 infection and inflammation. The objective was to provide a quantitative explanation for the range of COVID-19 symptom severity among individuals and to reveal the discriminatory importance of modeled mechanisms. The hypothesis we explored was that high levels of inflammation in COVID-19 may produce a significant amount of damage to uninfected cells. These cells would then produce DAMPs that further stimulate the inflammatory response. This model produced predictions that are consistent with clinical observations.
Hyperinflammation and disease severity. The model predicts that those with higher Disease Scores had substantially higher levels of inflammation (Fig. 2a). Further, peak inflammation was not reached until much later in those simulations leading to hyperinflammation (Fig. 3h). These modeling results are consistent with clinical findings that non-survivors tend to experience hyperinflammation and increasing levels of inflammatory biomarkers up to the time of death 7,9,47 . Our analysis provides an explanation of these findings.
In particular, the numerical bifurcation analysis suggests that viral infection can push the immune system into a self-sustaining high inflammatory state that can persist well beyond the resolution of infection. This hyperinflammatory state causes additional damage to uninfected cells, consequently leading to much higher Disease Scores. Due to the strong association between hyperinflammation and the Disease Score, we used those simulations that led to hyperinflammation as an in silico description of severe COVID-19. The percentage of simulations exhibiting hyperinflammation (13.5%) closely matched the reported proportion of the infected population with severe COVID-19 (14%), further motivating this decision 1 .

Viral load dynamics and disease severity.
In terms of viral load, the main difference that was reported between severe COVID-19 and those with milder disease pertained to the slope of VL after disease onset 7 . In particular, a faster VL decay was observed for those with milder disease after symptom onset, as measured using nasopharyngeal swabs 7 . This finding was also reported for VL sampled using other means and from various physiological compartments 48 . Interestingly, disease severity does not correlate strongly with peak viral loads 7,48 .
Our simulation results are consistent with these findings (Fig. 3b and c).The model offers an explanation to the slower viral load decay after peak infection in cases of severe COVID-19 through parameter κ 1 ; a lower κ 1 value both results in slower clearance of productively infected cells (see Eq. 3) and increases the risk of hyperinflammation (see Fig. 5). In fact, having a low κ 1 value was the most important predictor of hyperinflammation (Fig. 5b). In the model, low κ 1 values lead to the persistence of I and J cells, thereby prolonging PAMP and DAMP signaling and its downstream impact on the inflammatory response. The heightened inflammatory response promotes the generation of damaged cells J , strengthening DAMP signaling. Our analysis suggests that the ensuing feedback loop is the hallmark of hyperinflammation and severe COVID-19.
Characteristics of the innate immune response and disease severity. Many clinically observed associations were found between components of the innate immune response and disease severity. Among them, a dysregulated IFN response has been repeatedly associated with severe COVID-19 10,49 . Comparatively, our model suggests that a weaker IFN antiviral response (lower κ 0 ) leads to a higher likelihood of hyperinflammation (see Figs. 4c, 5a,b). One of the roles of type-I IFN is to limit the number of target cells that can be infected 31 .
Although it did not have a big impact on peak VL or the time to reach this peak in our simulations, an efficient type-I IFN response (higher κ 0 values) had a significant impact on inflammation. IFN by limiting the rate and the shear number of cells that carry PAMPs also constrains inflammation. Poor NK-cell cytotoxic ability was also linked with severe COVID-19 50 . NK-cells are important actors in the innate immune response that can target both infected and damaged cells and release molecules that precipitate their apoptosis 29,30 . In the model, poorer clearance of cells having PAMPs and DAMPs by the innate immune response is represented by lower κ 1 values, the most important predictor of hyperinflammation and severe COVID-19 (Fig. 5b). A poorer response from NK cells could hence lead to severe COVID-19 by enabling prolonged PAMP and DAMP signaling. Patients with severe COVID-19 also have increased levels of neutrophil extracellular traps (NETs) 25,26 . These can lead to immunothrombosis and the generation of damaged cells 23,24 . In the model, higher generation of damaged cells is represented by larger values of the parameter ν (Eq. 5). The bifurcation analysis suggests that larger values of ν lowers the threshold for the number of damaged cells that are required to reach a hyperinflammatory state (Fig. S1), facilitating severe COVID-19.
Finally, clinically defined severe COVID-19 has been associated with higher abundance and activation of proinflammatory macrophages 51 . In the model, those with higher disease scores had a higher number of innate immune cells ( D 0 + D 1 ) and a higher proportion of these cells were activated. The parameters that dictate innate immune cell recruitment ( ) and activation ( σ ) were both strongly associated with hyperinflammation and having a high Disease Score. www.nature.com/scientificreports/ clinical findings, suggesting that the model we developed here will be a useful tool to understand SARS-CoV-2 pathogenesis and predicting the impact of treatment. Indeed, observational studies revealed the use of corticosteroids can lead to more severe symptoms in those with milder disease, but generally improves outcomes for those with more severe symptoms [51][52][53][54][55] . Further, it is reported that those with milder COVID-19 generally experience slower viral load decay under corticosteroid treatment, an effect that was not found to be statistically significant among those with severe disease 52 . This latter result was also observed in our simulations (Fig. S8). One of the roles of inflammation is to stimulate cytolytic and phagocytic activities. By lowering this ability among those who experienced milder disease, the use of corticosteroids may lead to slower viral clearance. Hence, corticosteroids could have both a negative effect (slower clearance of PAMP carrying cells) and a beneficial effect (slower rate of bystander cell damage) on disease pathogenesis in those with mild disease. This could result in corticosteroids sometimes improving, sometimes worsening disease severity. Comparatively, the net beneficial effect of corticosteroids on those with severe COVID-19 may be the result of a smaller downstream impact of the drug on already less efficient NK cells 50 . For antivirals, simulations suggest that early administration is crucial for antivirals to impact disease severity, particularly for those that would have experienced severe COVID-19 in absence of treatment (Fig. 6). Our results also suggest the existence of an inflammation threshold beyond which antivirals may be unable to prevent hyperinflammation. It also suggests that early administration could reduce disease severity by preventing the infection from driving inflammation across the threshold. Comparatively, many clinical trials failed to show a substantial effect of antivirals 56 . More recently, a reduction of around 50% in the risk of hospitalization was observed after administering the antiviral molnupiravir 57 . The studied cohort consisted of individuals that had mild/moderate symptoms and were not expected to be hospitalized within 48 h of randomization 57 . Our results suggest this latter criteria may be crucial to ensure the beneficial effect of antivirals on disease severity.

Effect of corticosteroids and antivirals
Fitting clinical markers of disease severity. We used the disease score as a proxy measure of disease severity in our model. Quantitatively fitting our model by comparing the disease score with usual clinical disease severity scores may be inappropriate as it would require an established relationship between these quantities. In addition, clinical score data may be too sparse to reliably estimate the parameters in our model. Instead, biologically plausible parameter ranges were used to investigate gross changes in disease severity. Accordingly, changes that give rise to sufficiently different distributions of disease severity, such as those effected by some variants of concern, may inform parameters characterizing the variants. For example, the delta variant may cause overall higher viral loads 58 , which would lead to stronger immune activation on average, explaining the higher severity of this variant.
Other within-host models of SARS-CoV-2 infections. The model we propose is unique, but many of the effects it includes are present in other within-host models of SARS-CoV-2 infection. Despite differences in model formulations, there is a general agreement across models about the necessity for early antiviral administration [12][13][14][15] . However, models differed in terms of the components of the inflammatory, innate or adaptive immune response they include [11][12][13][14][15][16][17] . Some of the models included an effect of type-I interferon, either as a promoter of the pro-inflammatory response or helping susceptible cells resist infection [11][12][13] . One of these models concluded, as per our analysis, that an inefficient type-I IFN response can lead to accentuated tissue damage 11 . In comparison, our model simultaneously explains the more important clinical associations between severe disease and biomarkers of the innate immune response. It distinguishes itself by the inclusion of both a positive feedback loop between damaged cells and inflammation and an effect of innate immune cells on both infected and damaged cells. This latter effect highlights the importance of cells capable of clearing both types of cells in the pathogenesis of severe COVID-19.
Clinically, it has been shown that many individual characteristics, such as age and immune markers, e.g. lymphocyte counts and neutrophil counts, are associated with COVID-19 disease severity 59 .This complexity contrasts with the simplicity of the model we propose. While it is likely that many factors are involved in severe disease development, our model represents a high-level description of some of the physiological pathways, in particular, the innate immune response pathways, that may be involved in disease progression. This is consistent with the data from Li et al. and Lasso et al. that many immune markers predict disease severity better than other features, such as demographic characteristics 59,60 . Although we investigated more complex models, we decided against the modeling of individual cytokines or cells, as they often exhibit overlapping functions, and because some of these functions have been poorly studied leading to uncertainties in parameter values. However, some of the more complex models reported in the literature did give rise to interesting hypotheses that may warrant further investigation, such as the role of monocyte-to-macrophage differentiation or the role of anti-inflammatory cytokines 11,12 . Further, more complex models have provided complementary insight on the conditional impact of immunomodulatory treatment on disease severity that are consistent with the results reported herein 61,62 .

Conclusion
Our analysis revealed key aspects of the innate immune response that dictate inflammation trajectories and disease severity. The most important parameter suggested by bifurcation and decision tree analyses was κ 1 , representing the ability of the system to rid itself of cells carrying PAMPs and DAMPs. The analysis suggested that when this parameter is high enough, hyperinflammation can be avoided. The bifurcation analysis also suggested that small values of ν , representing the amount of bystander cell damage due to inflammation has a similar effect. In other words, the ability of the innate immune response to target PAMPs and DAMPs carrying cells appears to be key to the determination of COVID-19 severity. This suggests that therapies that specifically target aspects of the innate immune response may prove beneficial in comparison to broadly acting anti-inflammatory agents. www.nature.com/scientificreports/ The model also underlined the role of DAMPs in maintaining high levels of inflammation later in the course of infection. It suggests DAMPs may be an interesting therapeutic target for COVID-19. When exploring such novel treatment strategies, the model presented here could provide a means of exploring timing of treatment and dose effects in silico. Hopefully, a better understanding of the pathology of SARS-CoV-2 will lead to decreased mortality in this and similar diseases.